home *** CD-ROM | disk | FTP | other *** search
/ Pascal Super Library / Pascal Super Library (CW International)(1997).bin / MATH / NRPAS13 / SVDFIT.PAS < prev    next >
Pascal/Delphi Source File  |  1991-04-29  |  2KB  |  49 lines
  1. PROCEDURE svdfit(x,y,sig: glndata; ndata: integer; VAR a: glmma; mma: integer;
  2.       VAR u: glmpbynp; VAR v: glnpbynp; VAR w: glnparray; mp,np:
  3.       integer; VAR chisq: real);
  4. (* Programs using routine SVDFIT must define the types
  5. TYPE
  6.    glndata = ARRAY [1..ndata] OF real;
  7.    glmma = ARRAY [1..mma] OF real;
  8.    glmpbynp = ARRAY [1..mp,1..np] OF real;
  9.    glnpbynp = ARRAY [1..np,1..np] OF real;
  10.    glnparray = ARRAY [1..np] OF real;
  11. in the main routine. Implementations without conformant arrays require mma=np
  12. and the declaration
  13.    glnparray = glmma;
  14. and also mp=ndata, with the declaration
  15.    glmparray = glndata;
  16. in the main routine for use by the procedure SVBKSB. All user programs must also
  17. include a routine that computes the desired basis functions with the declaration
  18. PROCEDURE func(x: real: VAR p: glmma; mma: integer);
  19. which should return the values of first mma basis functions at x in p.
  20. FPOLY and FLEG (renamed to FUNC) are possible choices. *)
  21. CONST
  22.    tol=1.0e-5;
  23. VAR
  24.    j,i: integer;
  25.    wmax,tmp,thresh,sum: real;
  26.    b: glndata;
  27.    afunc: glmma;
  28. BEGIN
  29.    FOR i := 1 TO ndata DO BEGIN
  30.       func(x[i],afunc,mma);
  31.       tmp := 1.0/sig[i];
  32.       FOR j := 1 TO mma DO u[i,j] := afunc[j]*tmp;
  33.       b[i] := y[i]*tmp
  34.    END;
  35.    svdcmp(u,ndata,mma,mp,np,w,v);
  36.    wmax := 0.0;
  37.    FOR j := 1 TO mma DO IF (w[j] > wmax) THEN wmax := w[j];
  38.    thresh := tol*wmax;
  39.    FOR j := 1 TO mma DO IF (w[j] < thresh) THEN w[j] := 0.0;
  40.    svbksb(u,w,v,ndata,mma,mp,np,b,a);
  41.    chisq := 0.0;
  42.    FOR i := 1 TO ndata DO BEGIN
  43.       func(x[i],afunc,mma);
  44.       sum := 0.0;
  45.       FOR j := 1 TO mma DO sum := sum+a[j]*afunc[j];
  46.       chisq := chisq+sqr((y[i]-sum)/sig[i])
  47.    END
  48. END;
  49.